Data collection
### Pollution
Stations <- read_csv(here::here("ass_data",
"csv",
"Stacje.csv"))
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Nr = col_double(),
## `Kod stacji` = col_character(),
## `Kod międzynarodowy` = col_character(),
## `Nazwa stacji` = col_character(),
## `Stary Kod stacji` = col_character(),
## `Data uruchomienia` = col_character(),
## `Data zamknięcia` = col_character(),
## `Typ stacji` = col_character(),
## `Typ obszaru` = col_character(),
## `Rodzaj stacji` = col_character(),
## Województwo = col_character(),
## Miejscowość = col_character(),
## Ulica = col_character(),
## `WGS84 φ N` = col_double(),
## `WGS84 λ E` = col_double()
## )
# Following piece of code is commented out because running it might take up to 1h
# The output of this section of the code is the data frame called "pollution_yearly"
# which consist annual reading of the eight pollutants ("pm10","pm25","so2","o3","nox","no2","co","c6h6")
# from the period of 2000-2019 for all the 182 subregions in malopolska
#### Loading all pollutants 2000-2019 data
pollutants <- c("pm10","pm25","so2","o3","nox","no2","co","c6h6")
pollutants_csv <- list()
for(i in pollutants){
#loading csv for whole country and slicing first row
filepath <- file.path(here::here("ass_data","csv"),paste(toupper(i),".csv",sep=""))
assign(sprintf("%s_poland",i), slice(read_csv(filepath),-1))
pollutants_csv[[match(i,pollutants)]] <- sprintf("%s_poland",i)
}
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Rok = col_character(),
## Województwo = col_character(),
## `Kod strefy` = col_character(),
## `Nazwa strefy` = col_character(),
## `Kod stacji` = col_character(),
## Wskaźnik = col_character(),
## `Czas uśredniania` = col_character(),
## Średnia = col_character(),
## Min = col_character(),
## Maks = col_character(),
## `L>50 (S24)` = col_character(),
## `36 maks. (S24)` = col_character(),
## `Perc. 90.4 (S24)` = col_character(),
## `Maks (S24)` = col_character(),
## `Liczba pomiarów` = col_character(),
## Kompletność = col_character(),
## `Liczba Lato/Zima` = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Rok = col_character(),
## Województwo = col_character(),
## `Kod strefy` = col_character(),
## `Nazwa strefy` = col_character(),
## `Kod stacji` = col_character(),
## Wskaźnik = col_character(),
## `Czas uśredniania` = col_character(),
## Średnia = col_character(),
## Min = col_character(),
## Maks = col_character(),
## `Liczba pomiarów` = col_character(),
## Kompletność = col_character(),
## `Liczba Lato/Zima` = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character()
## )
## ℹ Use `spec()` for the full column specifications.
##
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character()
## )
## ℹ Use `spec()` for the full column specifications.
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Rok = col_character(),
## Województwo = col_character(),
## `Kod strefy` = col_character(),
## `Nazwa strefy` = col_character(),
## `Kod stacji` = col_character(),
## Wskaźnik = col_character(),
## `Czas uśredniania` = col_character(),
## Średnia = col_character(),
## Min = col_character(),
## Maks = col_character(),
## `Liczba pomiarów` = col_character(),
## Kompletność = col_character(),
## `Liczba Lato/Zima` = col_character()
## )
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Rok = col_character(),
## Województwo = col_character(),
## `Kod strefy` = col_character(),
## `Nazwa strefy` = col_character(),
## `Kod stacji` = col_character(),
## Wskaźnik = col_character(),
## `Czas uśredniania` = col_character(),
## Średnia = col_character(),
## Min = col_character(),
## Maks = col_character(),
## `L>200 (S1)` = col_character(),
## `19 maks. (S1)` = col_character(),
## `Perc. 99.8 (S1)` = col_character(),
## `Liczba pomiarów` = col_character(),
## Kompletność = col_character(),
## `Liczba Lato/Zima` = col_character()
## )
## Warning: Missing column names filled in: 'X15' [15], 'X16' [16], 'X17' [17],
## 'X18' [18], 'X19' [19], 'X20' [20]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## X15 = col_logical(),
## X16 = col_logical(),
## X17 = col_logical(),
## X18 = col_logical(),
## X19 = col_logical(),
## X20 = col_logical()
## )
## ℹ Use `spec()` for the full column specifications.
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Rok = col_character(),
## Województwo = col_character(),
## `Kod strefy` = col_character(),
## `Nazwa strefy` = col_character(),
## `Kod stacji` = col_character(),
## Wskaźnik = col_character(),
## `Czas uśredniania` = col_character(),
## Średnia = col_character(),
## Min = col_character(),
## Maks = col_character(),
## `Liczba pomiarów` = col_character(),
## Kompletność = col_character(),
## `Liczba Lato/Zima` = col_character()
## )
# ##### Joining pollution readings with the station data and converting it to SF
#
# for (i in pollutants_csv) {
# # joining with stations dataframe by station code
# filename <- eval(parse(text = i))
# assign(sprintf("%s_sf",i), left_join(filename,Stations,
# by = c("Kod stacji" = "Kod stacji")))
# # droping na from "WGS84 φ N" column
# filename_2 <- eval(parse(text = sprintf("%s_sf",i)))
# assign(sprintf("%s_sf",i), drop_na(filename_2,"WGS84 φ N") )
# # converting to SF
# filename_2 <- eval(parse(text = sprintf("%s_sf",i)))
# assign(sprintf("%s_sf",i), st_as_sf(filename_2, coords = c("WGS84 λ E","WGS84 φ N"), crs = 4326))
# }
#
# ##### putting all pollutants SFs to a list
# pollutants_sf<- list(pm10_poland_sf,pm25_poland_sf,co_poland_sf,no2_poland_sf,
# nox_poland_sf,so2_poland_sf,c6h6_poland_sf,o3_poland_sf)
# #### Interpolation
#
# ##### SP outline od Malopolska
# ls_border <- lesserpoland %>%
# st_transform(., 2180) %>%
# as(., 'Spatial')
#
# ##### SF -> SP yearly measurments - All pollutants
# pollutants_yearly <- lapply(pollutants_sf,function(x) group_split(x,Rok) )
#
# pollutants_yearly_sp <- lapply(pollutants_yearly,
# function(x) lapply(x,
# function(y) as(st_transform(y, 2180), 'Spatial')))
#
# ##### creating a grind
# emptygrd <- as.data.frame(spsample(ls_border, n=1000, type="regular", cellsize=1000))
#
# names(emptygrd) <- c("X", "Y")
#
# coordinates(emptygrd) <- c("X", "Y")
#
# gridded(emptygrd) <- TRUE # Create SpatialPixel object
# fullgrid(emptygrd) <- TRUE # Create SpatialGrid object
#
# ##### creating raster with yearly measurments from 2000-2019 for all pollutants
#
# pollutants_yearly_ras<- list()
#
# for (p in pollutants_yearly_sp) {
# # Add the projection to the grid
# proj4string(emptygrd)<- proj4string(p[[1]])
# # Interpolate the grid cells using a power value of 2
# assign("interpolate_list", lapply(p,function(i) gstat::idw(Średnia ~ 1, i, newdata=emptygrd, idp=2.0)))
# # Convert output to raster object
# assign("ras", lapply(interpolate_list,function(i) raster(i)))
# # Clip the raster to Lesserpoland outline
# assign("ras_masked" ,lapply(ras,function(i) raster::mask(i, ls_border)))
# # add to the pollutants_yearly_ras
# pollutants_yearly_ras<-append(pollutants_yearly_ras,list(ras_masked))
# }
#
# ##### Plot the raster - test
# # plot(pollutants_yearly_ras[[8]][[20]])
#
# ##### extracting values FOR EVERY REGION from a raster and adding it to dataframe
# pollutants_yearly_mean<- list()
#
# for (p in pollutants_yearly_ras) {
# assign("pollutants_yearly_gminy", lapply(p, function(y)
# lapply(gminy_LP_grouped, function(g)
# raster::extract(y,g,small=TRUE))))
# pollutants_yearly_mean<-append(pollutants_yearly_mean,list(pollutants_yearly_gminy))
# }
#
# ##### calculation annual mean for every region
# pollutants_yearly_gminy_mean<- pollutants_yearly_mean %>%
# lapply(., function(p)
# lapply(p,function(y)
# lapply(y,function(g)
# mean(as.numeric(as.character(unlist(g))))
# )))
#
# #####transfering to data frame
#
# gminy_codes <- gminy_LP %>%
# dplyr::select(.,"JPT_KOD_JE","JPT_NAZWA_") %>%
# st_drop_geometry() %>%
# mutate(JPT_KOD_JE = as.numeric(JPT_KOD_JE)) %>%
# arrange(.,JPT_KOD_JE)
#
# gminy_codes_v<- as.vector(gminy_codes$JPT_KOD_JE)
#
# ##### naming list elements + convert to the final data frame
# df <- pollutants_yearly_gminy_mean %>%
# lapply(., function(p)
# lapply(p,function(y)
# as.data.frame(y))) %>%
# lapply(., function(p)
# lapply(p, function(y) y %>%
# rename_at(vars(colnames(y)), function(x) gminy_codes_v) %>%
# pivot_longer(everything())))
#
# names(df) <- c("pm10","pm25","co","no2","nox","so2","c6h6","o3")
#
# for (i in c(1,4,5,6,8)) {
# names(df[[i]]) <- as.character(unlist(seq(2000,2019,by=1)))
# names(df[[2]]) <- as.character(unlist(seq(2002,2019,by=1)))
# names(df[[7]]) <- as.character(unlist(seq(2001,2019,by=1)))
# names(df[[3]]) <- as.character(unlist(seq(2003,2019,by=1)))
# }
#
#
# df_pollutants <- df %>%
# lapply(., function(x) x %>% reduce(left_join, by = "name"))
#
# for (i in c(1,4,5,6,8)) {
# names(df_pollutants[[i]]) <- c("area_code",as.character(unlist(seq(2000,2019,by=1))))
# names(df_pollutants[[2]]) <- c("area_code",as.character(unlist(seq(2002,2019,by=1))))
# names(df_pollutants[[7]]) <- c("area_code",as.character(unlist(seq(2001,2019,by=1))))
# names(df_pollutants[[3]]) <- c("area_code",as.character(unlist(seq(2003,2019,by=1))))
# }
#
# df_pollutants <- df_pollutants %>%
# reduce(left_join, by = "area_code",all.x = TRUE,all.y = TRUE)
#
#
# df_pollutants <- df_pollutants%>%
# rename_at(.vars = vars(seq(136,155,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_o3", .))) %>%
# rename_at(.vars = vars(seq(117,135,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_c6h6", .))) %>%
# rename_at(.vars = vars(seq(97,116,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_so2", .)))%>%
# rename_at(.vars = vars(seq(77,96,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_nox", .))) %>%
# rename_at(.vars = vars(seq(57,76,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_no2", .)))%>%
# rename_at(.vars = vars(seq(40,56,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_co", .))) %>%
# rename_at(.vars = vars(seq(22,39,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_pm25", .)))%>%
# rename_at(.vars = vars(seq(2,21,1)),
# .funs = funs(sub("[.x]*$|[.y]*$", "_pm10", .)))
# ##### saving to csv
# write.csv(df_pollutants,here::here("ass_data",
# "csv",
# "pollution_yearly.csv"), row.names = FALSE)
#### Loading generated csv file with pollutans
pollutants <- read_csv(here::here("ass_data",
"csv",
"pollution_yearly.csv"))
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
#### calculating mean reading of the pollutants for the whole region of Malopolska
pollutants_ls_mean<- pollutants%>%
summarise_all( mean) %>%
pivot_longer(-area_code,
names_to = "year_pollutant",
values_to = "value",
) %>%
dplyr::select(-area_code) %>%
mutate(year= as.factor(str_extract(year_pollutant, "^\\d+")),
pollutant= as.character(str_extract(year_pollutant, "(?<=_)(.*)")))
#### calculating difference of the annual PM25 reading before 2008 and after
pm25 <- pollutants%>%
dplyr::select(c(1,22:39))%>%
mutate(across(.cols = 2:19, as.numeric)) %>%
mutate(pm25_08_18= round(rowMeans(.[c(9:18)], na.rm = TRUE),1)) %>%
mutate(pm25_02_08= round(rowMeans(.[c(2:8)], na.rm = TRUE),1)) %>%
mutate(pm25= pm25_08_18) %>%
mutate(pm25_dif= pm25_08_18-pm25_02_08) %>%
dplyr::select(area_code,pm25,pm25_08_18,pm25_02_08,pm25_dif)
### Density of population
density <- read_csv(here::here("ass_data",
"csv",
"density.csv")) %>%
slice(.,-2) %>%
row_to_names(row_number = 1) %>%
dplyr::rename(.,area_code=1,area_name=2) %>%
right_join(.,gminy_codes, by=c("area_code"="JPT_KOD_JE"))
## Warning: Missing column names filled in: 'X4' [4], 'X5' [5], 'X6' [6], 'X7' [7],
## 'X8' [8], 'X9' [9], 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14],
## 'X15' [15], 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## Kod = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## Warning in row_to_names(., row_number = 1): Row 1 does not provide unique names.
## Consider running clean_names() after row_to_names().
#### calculating difference between 02-08 and 08-19
density_08_19 <- density %>%
dplyr::select(c(1,3:20))%>%
mutate(across(.cols = 2:19, as.numeric)) %>%
mutate(density_km2_08_19= round(rowMeans(.[c(9:19)], na.rm = TRUE),0)) %>%
mutate(density_km2= density_km2_08_19) %>%
dplyr::select(area_code,density_km2)
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 2:19, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
### population
pop <- read_csv(here::here("ass_data",
"csv",
"ludnosc.csv")) %>%
slice(.,c(-1,-2,-3,-5)) %>%
row_to_names(row_number = 1) %>%
dplyr::rename(.,area_code=1,area_name=2) %>%
mutate(across(.cols = 3:22, as.numeric))%>%
mutate(population_00_08 = round(rowMeans(subset(., select = as.numeric(c(3:11))), na.rm = TRUE),0)) %>%
mutate(population_08_19 = round(rowMeans(subset(., select = as.numeric(c(12:22))), na.rm = TRUE),0)) %>%
mutate(population_dif = population_08_19-population_00_08 )
## Warning: Missing column names filled in: 'X4' [4], 'X5' [5], 'X6' [6], 'X7' [7],
## 'X8' [8], 'X9' [9], 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14],
## 'X15' [15], 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20],
## 'X21' [21], 'X22' [22]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## Kod = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## Warning in row_to_names(., row_number = 1): Row 1 does not provide unique names.
## Consider running clean_names() after row_to_names().
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
pop_dif <- pop %>%
dplyr::select(area_code,population_dif) %>%
mutate(population_dif = replace_na(population_dif, 0))
pop_08_19 <- pop %>%
dplyr::select(area_code,population_08_19)
### Furnaces removed in malopolska during 2008-2019
furnaces <- read_csv(here::here("ass_data",
"csv",
"furnaces_final.csv")) %>%
dplyr::select(area_code,area_name,podregion,powiat_name,
typ_gminy,'2008','2009','2010','2011','2012','2013','2014','2015','2016','2017','2018') %>%
drop_na(area_code)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## area_code = col_double(),
## area_name = col_character(),
## podregion = col_character(),
## powiat_name = col_character(),
## typ_gminy = col_character(),
## `2008` = col_double(),
## `2009` = col_double(),
## `2010` = col_double(),
## `2011` = col_double(),
## `2012` = col_double(),
## `2013` = col_double(),
## `2014` = col_double(),
## `2015` = col_double(),
## `2016` = col_double(),
## `2017` = col_double(),
## `2018` = col_double()
## )
furnaces_final <- furnaces %>%
mutate(type_n=case_when(str_detect(typ_gminy, "City|Urban$") ~ "1",
str_detect(typ_gminy, "Urban-rural") ~ "3",
str_detect(typ_gminy, "Rural") ~ "2")) %>%
mutate(area_code=paste(area_code, type_n,sep="")) %>%
mutate(all_removed= rowSums(.[sapply(.,is.numeric)], na.rm = TRUE)) %>%
mutate(mean_removed= round(rowMeans(.[sapply(.,is.numeric)], na.rm = TRUE),0)) %>%
mutate(area_code=as.numeric(area_code)) %>%
left_join(.,pop_08_19,by="area_code") %>%
mutate(furnaces_per_1000= round(all_removed/(population_08_19/1000),0))
total_fur <- sum(furnaces_final$all_removed)
furnaces_08_19 <- furnaces_final %>%
dplyr::select(area_code,furnaces_per_1000)
gminy_LP_fur <- gminy_LP %>%
mutate(JPT_KOD_JE=as.numeric(JPT_KOD_JE)) %>%
left_join(.,furnaces_final,
by = c("JPT_KOD_JE"="area_code"))
### roads length in subregions
# following piece of code is commented out because running it might take up to 20min
# The output of this section of the code is the data frame called "roads_length"
# which consist of total road length in metres for all the 182 subregions in malopolska
# calculated using OSM data
# roads <- st_read(here::here("ass_data",
# "geo", 'malopolskie-latest-free',
# "gis_osm_roads_free_1.shp"))%>%
# st_transform(., 2180)
#
# roads_length <- lapply(gminy_LP_grouped,function(x) sum(st_length(sf::st_intersection(roads,x))))
#
# names<- gminy_codes%>%
# arrange(., JPT_KOD_JE) %>%
# dplyr::select(JPT_KOD_JE) %>%
# pull()
#
#
# roads_df<-as.data.frame(roads_length)
# colnames(roads_df)<- names
#
# roads_df<- roads_df %>%
# pivot_longer(everything(),
# names_to="area_code",
# values_to="roads_length_m") %>%
# mutate(roads_length_m = as.numeric(gsub("[m]$", "", roads_length_m))) %>%
# mutate(area_code=as.numeric(area_code))%>%
# left_join(.,gminy_codes,by=c("area_code"="JPT_KOD_JE"))
#
# write.csv(roads_df,here::here("ass_data",
# "csv",
# "roads_length.csv"), row.names = FALSE)
roads <- read_csv(here::here("ass_data",
"csv",
"roads_length.csv")) %>%
mutate(roads_km=round(roads_lenght_m/1000,3))
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## area_code = col_double(),
## roads_lenght_m = col_double(),
## JPT_NAZWA_ = col_character()
## )
### forrest area
forrest <- read_csv(here::here("ass_data",
"csv",
"forrest_area_prc.csv")) %>%
slice(.,c(-2))%>%
row_to_names(row_number = 1) %>%
dplyr::rename(.,area_code=1,area_name=2) %>%
mutate(across(.cols = 3:22, as.numeric)) %>%
left_join(gminy_codes,.,by=c("JPT_KOD_JE"="area_code")) %>%
mutate(forrest_08_19 = round(rowMeans(subset(., select = as.numeric(c(13:23))), na.rm = TRUE),0)) %>%
mutate(forrest = forrest_08_19) %>%
dplyr::select(JPT_KOD_JE,forrest)
## Warning: Missing column names filled in: 'X4' [4], 'X5' [5], 'X6' [6], 'X7' [7],
## 'X8' [8], 'X9' [9], 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14],
## 'X15' [15], 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20],
## 'X21' [21], 'X22' [22]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## Kod = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## Warning in row_to_names(., row_number = 1): Row 1 does not provide unique names.
## Consider running clean_names() after row_to_names().
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:22, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
### land use data
# following piece of code is commented out because running it might take up to 20min
# The output of this section of the code is the data frame called "urban_area"
# which consist of percentages of urban development in the subregion for all the 182 subregions in malopolska
# calculated using OSM data
# land <- st_read(here::here("ass_data",
# "geo", 'malopolskie-latest-free',
# "gis_osm_landuse_a_free_1.shp"))%>%
# st_transform(., 2180)
#
# urban<-land %>%
# dplyr::filter(fclass=="industrial"|fclass=="retail"|fclass=="residential"|fclass=="commercial"|fclass=="heath")
#
# urban_area <- lapply(gminy_LP_grouped,function(x) sum(st_area(sf::st_intersection(urban,x))))
# total_area <- lapply(gminy_LP_grouped,function(x) st_area(x))
#
#
# urban_df<-as.data.frame(urban_area)
# total_df<- as.data.frame(total_area)
#
# colnames(urban_df)<- names
# colnames(total_df)<- names
#
# urban_df<- urban_df %>%
# pivot_longer(everything(),
# names_to="area_code",
# values_to="urban_area_m2") %>% left_join(.,
# pivot_longer(total_df,everything(),
# names_to="area_code",
# values_to="total_area_m2"),
# by="area_code") %>%
# mutate(total_area_m2=as.numeric(total_area_m2),
# urban_area_m2=as.numeric(urban_area_m2),
# area_code=as.numeric(area_code))%>%
# left_join(.,
# gminy_codes,
# by=c("area_code"="JPT_KOD_JE")) %>%
# mutate(urban_landuse = urban_area_m2/total_area_m2*100) %>%
# mutate_at(vars(-area_code,-JPT_NAZWA_), funs(round(., 1)))
#
# write.csv(urban_df,here::here("ass_data",
# "csv",
# "urban_area.csv"), row.names = FALSE)
urban<-read_csv(here::here("ass_data",
"csv",
"urban_area.csv")) %>%
dplyr::select(area_code,urban_landuse,total_area_m2)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## area_code = col_double(),
## urban_area_m2 = col_double(),
## total_area_m2 = col_double(),
## JPT_NAZWA_ = col_character(),
## urban_landuse = col_double()
## )
### altitude data
# following piece of code is commented out because running it might take up to 20min
# The output of this section of the code is the data frame called "elevation."
# which consist of average altitude for all the 182 subregions in malopolska
# calculated using elevatr package
# ls_border <- lesserpoland %>%
# st_transform(., 2180) %>%
# as(., 'Spatial')
#
# #### creating the dataframe with SP bbox
# loc_df <- data.frame(x = runif(6,min=sp::bbox(buffer(ls_border,10000))[1,1],
# max=sp::bbox(buffer(ls_border,10000))[1,2]),
# y = runif(6,min=sp::bbox(buffer(ls_border,10000))[2,1],
# max=sp::bbox(buffer(ls_border,10000))[2,2]))
#
# #### getting the altitude data using the elevatr package
# x <- get_elev_raster(locations = loc_df, prj = sp::proj4string(ls_border), z=10)
#
# elev_mask<-mask(x,ls_border)
#
# writeRaster(elev_mask, filename=here::here("ass_data",
# "csv",
# "ls_elevation.tif"), format="GTiff", overwrite=TRUE)
#
# #### extracting mean altitude from a raster and adding it to dataframe
#
# elev_gminy <- lapply(gminy_LP_grouped,function(x) raster::extract(elev_mask,x,small=TRUE))
# elev_gminy_mean <- elev_gminy %>%
# lapply(.,function(x) as.numeric(as.character(unlist(x)))) %>%
# lapply(.,function(x) mean(x,na.rm = TRUE))
#
# plot(mask(elev_mask,gminy_LP_grouped[[180]]))
#
# elev_df<-as.data.frame(elev_gminy_mean)
#
# colnames(elev_df)<- names
#
# elev_df<- elev_df %>%
# pivot_longer(everything(),
# names_to="area_code",
# values_to="elevation_m") %>%
# mutate(area_code=as.numeric(area_code))%>%
# left_join(.,
# gminy_codes,
# by=c("area_code"="JPT_KOD_JE")) %>%
#
# write.csv(elev_df,here::here("ass_data",
# "csv",
# "elevation.csv"), row.names = FALSE)
#### loading the altitude data
altitude<-read_csv(here::here("ass_data",
"csv",
"elevation.csv")) %>%
rename(altitude_m_masl=elevation_m)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## area_code = col_double(),
## elevation_m = col_double(),
## JPT_NAZWA_ = col_character()
## )
#### gminy budget
budget<-read_csv(here::here("ass_data",
"csv",
"gminy_budget.csv")) %>%
slice(.,c(-1,-3))%>%
row_to_names(row_number = 1) %>%
dplyr::rename(.,area_code=1,area_name=2) %>%
mutate(across(.cols = 3:21, as.numeric)) %>%
left_join(gminy_codes,.,by=c("JPT_KOD_JE"="area_code")) %>%
mutate(budget_08_19 = round(rowMeans(subset(., select = as.numeric(c(11:22))), na.rm = TRUE),0)) %>%
mutate(budget = budget_08_19) %>%
dplyr::select(JPT_KOD_JE,budget)
## Warning: Missing column names filled in: 'X4' [4], 'X5' [5], 'X6' [6], 'X7' [7],
## 'X8' [8], 'X9' [9], 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14],
## 'X15' [15], 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20],
## 'X21' [21]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## Kod = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## Warning in row_to_names(., row_number = 1): Row 1 does not provide unique names.
## Consider running clean_names() after row_to_names().
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:21, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
#### gminy flats
flats<-read_csv(here::here("ass_data",
"csv",
"flats.csv")) %>%
slice(.,c(-2))%>%
row_to_names(row_number = 1) %>%
dplyr::rename(.,area_code=1,area_name=2) %>%
mutate(across(.cols = 3:20, as.numeric)) %>%
left_join(gminy_codes,.,by=c("JPT_KOD_JE"="area_code")) %>%
mutate(flats_08_19 = round(rowMeans(subset(., select = as.numeric(c(10:21))), na.rm = TRUE),1)) %>%
mutate(flats = flats_08_19) %>%
dplyr::select(JPT_KOD_JE,flats)
## Warning: Missing column names filled in: 'X4' [4], 'X5' [5], 'X6' [6], 'X7' [7],
## 'X8' [8], 'X9' [9], 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14],
## 'X15' [15], 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## Kod = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## Warning in row_to_names(., row_number = 1): Row 1 does not provide unique names.
## Consider running clean_names() after row_to_names().
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(.cols = 3:20, as.numeric)`.
## Warning in fn(col, ...): NAs introduced by coercion
Descriptive statistics
### assembling the final data frame used for the analysis
df <- gminy_LP %>%
dplyr::select(JPT_KOD_JE,JPT_NAZWA_) %>%
rename(area_code=JPT_KOD_JE,
area_name=JPT_NAZWA_) %>%
mutate(area_code= as.numeric(area_code)) %>%
left_join(.,density_08_19, by="area_code") %>%
left_join(.,furnaces_08_19,by="area_code") %>%
left_join(.,pm25, by="area_code") %>%
left_join(.,pop_08_19, by="area_code") %>%
left_join(.,roads[c("area_code","roads_km")], by="area_code") %>%
left_join(.,urban, by="area_code") %>%
left_join(.,altitude,by="area_code") %>%
left_join(.,forrest, by=c("area_code"="JPT_KOD_JE")) %>%
left_join(.,budget, by=c("area_code"="JPT_KOD_JE")) %>%
left_join(.,flats, by=c("area_code"="JPT_KOD_JE")) %>%
mutate(roads_density_km2= roads_km/(total_area_m2/1000000))
write<- df %>%
st_drop_geometry() %>%
write.csv(.,here::here("ass_data",
"csv",
"df.csv"), row.names = FALSE)
df_tidy<- df %>%
st_drop_geometry() %>%
dplyr::select(is.numeric,-area_code)%>%
pivot_longer(everything(),
names_to="All_variables",
values_to="val")%>%
mutate(All_variables = tolower(All_variables))
## Warning: Predicate functions must be wrapped in `where()`.
##
## # Bad
## data %>% select(is.numeric)
##
## # Good
## data %>% select(where(is.numeric))
##
## ℹ Please update your code.
## This message is displayed once per session.
### plotting annual polluntion trends in the region of lesserpoland
trend_ls <- pollutants_ls_mean %>%
dplyr::filter(pollutant=="pm25") %>%
ggplot(., aes(x=year,y=value)) +
geom_line(aes(group=pollutant),size=.25)+
geom_point(size=2,colour="red")+ geom_vline(aes(xintercept="2008"),linetype = 2,size=.5)+
geom_text(aes(x=7.25, label="MAQP", y=21), colour="black", angle=90,size=3 )+
labs(title=expression("Annual conentrantion of PM"[2.5]*" in Małopolska"),subtitle = "",
x ="Years", y = expression("Annual PM"[2.5]* " Concentration [μg/m"^3*"]" ))+
theme_minimal()
trend_ls
#ggsave("trend_ls.png", trend_ls, width=9, height=4,dpi=600)
### plotting study area map
ls_stations<- Stations %>%
left_join(pm25_poland,.,by = c("Kod stacji" = "Kod stacji")) %>%
drop_na(.,"WGS84 φ N") %>%
st_as_sf(.,coords = c("WGS84 λ E","WGS84 φ N"), crs = 4326) %>%
st_transform(., 2180) %>%
st_intersection(.,lesserpoland)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ls <- gminy_LP %>%
mutate(area_code=as.numeric(JPT_KOD_JE)) %>%
left_join(.,density_08_19,by="area_code")
cities <- gminy_LP %>%
filter(str_detect(JPT_KOD_JE,'^126'))
cities$JPT_NAZWA_[cities$JPT_NAZWA_ == "Nowy Sącz"] <- "Nowy Sacz"
tmap_mode("plot")
## tmap mode set to plotting
ls_bb<- st_bbox(lesserpoland,
crs = st_crs(lesserpoland)) %>%
st_as_sfc()
main <- tm_shape(ls, bbbox = ls_bb) +
tm_polygons("density_km2",alpha=.75, palette="PuBu",style="fisher", title = 'Population Density\npeople per square km')+
tm_compass(type = "arrow", position = c("right", "bottom")) +
tm_shape(ls_stations)+ tm_dots(col="yellow",size = 0.05)+
tm_layout(legend.position = c("left","bottom"),
legend.text.size=0.6,
legend.title.size = 0.75,
legend.title.fontface = 1,
legend.height = 0.9,
frame=FALSE)+
tm_shape(cities)+ tm_borders(col = "black",lwd = 1.5)+
tm_symbols(col = "red", scale = .5)+
tm_text("JPT_NAZWA_", xmod=0, ymod=-1,size=0.7,fontface="bold",bg.color = "white")+
tm_layout(inner.margin=c(0.02,0.02,0.02,0.2))+
tm_credits("(a)", position=c(0,0.85), size=1.5)+
tm_add_legend('symbol',col="yellow",size = 0.35,
border.col = "black",
labels = "The air quality monitoring stations")
inset = tm_shape(wojewodztwa) + tm_fill(col="white")+ tm_borders(col = "black")+
tm_shape(lesserpoland) + tm_fill(col="lightblue")+
tm_borders(col = "black", lwd = 1)+
tm_layout(frame=FALSE,
bg.color = "transparent")+
tm_credits("(b)", position=c(0,0.85), size=1.5)
#tmap_save(main,insets_tm = inset,insets_vp=viewport(0.86, 0.7, width = 0.35, height = 0.35), filename="test.png", dpi=300)
### histogram dependent variable PM25
pm25_hist <- df_tidy%>%
filter(All_variables=="pm25")%>%
ggplot(., aes(x=val)) +
geom_histogram(aes(x = val, y = ..density..),color="black", fill="white")+
geom_density(colour="red", size=1, adjust=1)
pm25_hist
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
### Dependent variable descriptive statistics
summary(df$pm25)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.30 31.40 32.80 32.55 33.60 37.70
### histogram independent variable
independent_hist <- df_tidy%>%
filter(All_variables!="pm25")%>%
ggplot(., aes(x=val)) +
geom_histogram(aes(x = val, y = ..density..),color="black", fill="white")+
geom_density(colour="red", size=1, adjust=1)+
facet_wrap(~All_variables, scales = 'free')
independent_hist
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
### transforming data log1p independent
independent_log1phist <- df_tidy%>%
filter(All_variables!="pm25")%>%
ggplot(., aes(x=log1p(val))) +
geom_histogram(aes(x = log1p(val), y = ..density..),color="black", fill="white")+
geom_density(colour="red", size=1, adjust=1)+
facet_wrap(~All_variables, scales = 'free')
independent_log1phist
## Warning in log1p(val): NaNs produced
## Warning in log1p(val): NaNs produced
## Warning in log1p(val): NaNs produced
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 125 rows containing non-finite values (stat_bin).
## Warning: Removed 125 rows containing non-finite values (stat_density).
density_loghist <- df_tidy%>%
filter(All_variables=="density_km2")%>%
ggplot(., aes(x=log(val))) +
geom_histogram(aes(x = log(val), y = ..density..),color="black", fill="white")+
geom_density(colour="red", size=1, adjust=1)+
facet_wrap(~All_variables, scales = 'free')
density_loghist
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
### plotting histograms
supp.labs <- c("Log1p(Flats per 1000 residents)", "Forrest Area","PM2.5","Log1p(Roads Density)","Log1p(Altitude)")
names(supp.labs) <- c("flats", "forrest","pm25","roads_density_km2","altitude_m_masl")
Distribution_hist <-df %>%
st_drop_geometry() %>%
dplyr::select(is.numeric,-area_code)%>%
mutate(roads_density_km2=log1p(roads_density_km2),
altitude_m_masl=log1p(altitude_m_masl),
flats=log1p(flats)) %>%
pivot_longer(everything(),
names_to="All_variables",
values_to="val")%>%
mutate(All_variables = tolower(All_variables))%>%
filter(All_variables=="pm25"|
All_variables=="roads_density_km2"|
All_variables=="altitude_m_masl"|
All_variables=="flats"|
All_variables=="forrest")%>%
ggplot(., aes(x=val)) +
geom_histogram(aes(x = val, y = ..density..),color="black", fill="white")+
geom_density(colour="red", size=1, adjust=1)+
labs(x = NULL, y = NULL)+
facet_wrap(~All_variables,ncol=1,scales="free",dir="v" ,labeller =labeller(All_variables=supp.labs))+
theme_minimal()
Distribution_hist
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
### boxplots
pm25_box <- df_tidy%>%
filter(All_variables=="pm25")%>%
ggplot(., aes(y=val)) +
geom_boxplot()
pm25_box
log_furnaces_box <- df_tidy%>%
filter(All_variables=="furnaces_per_1000")%>%
ggplot(., aes(y=log1p(val))) +
geom_boxplot()
log_furnaces_box
### mapping pm25 and furnaces
tmap_mode("plot")
## tmap mode set to plotting
breaks_pm25 <- c(28,30,33.5,34,34.5,35,38)
pm25_08_18_map <- tm_shape(df, bbbox = ls_bb) +
tm_polygons(col='pm25_08_18',breaks=breaks_pm25,#style="jenks",
palette="YlOrRd",alpha=1,title=expression('PM'[2.5]* '[μg/m'^3*']'))+
tm_layout(title = "Period 2008-2018",
title.position = c('right','top'),
title.size = 1,
legend.position = c(0,0),
legend.text.size=0.6,
legend.title.size = 1,
legend.title.fontface = 1,
legend.height = 0.4,
frame=FALSE)+
tm_credits("(b)", position=c(0,0.85), size=1)
pm25_02_08_map <- tm_shape(df, bbbox = ls_bb) +
tm_polygons(col='pm25_02_08',breaks=breaks_pm25,#style="fisher",
palette="YlOrRd",alpha=1,title=expression('PM'[2.5]* '[μg/m'^3*']'))+
tm_layout(title = "Period 2002-2007",
title.position = c('right','top'),
title.size = 1,
legend.position = c(0,0),
legend.text.size=0.6,
legend.title.size = 1,
legend.title.fontface = 1,
legend.height = 0.4,
frame=FALSE)+
tm_credits("(a)", position=c(0,0.85), size=1)
pm25_dif_map <- tm_shape(df, bbbox = ls_bb) +
tm_polygons(col='pm25_dif',style="fisher",
palette="seq",alpha=1,title=expression('Change in PM'[2.5]))+
tm_layout(title = "Change (b)-(a)",
title.position = c('right','top'),
aes.palette = list(seq = "-Spectral"),
title.size = 1,
legend.position = c(0,0),
legend.text.size=0.6,
legend.title.size = 0.75,
legend.title.fontface = 1,
legend.height = 0.4,
frame=FALSE)+
tm_credits("(c)", position=c(0,0.85), size=1)
furnaces_map <- tm_shape(df, bbbox = ls_bb) +
tm_polygons(col='furnaces_per_1000',style="fisher",
palette="YlGnBu",alpha=1,title='Removed furnaces\nper 1000 residents')+
tm_layout(title = "Period 2008-2018",
title.position = c('right','top'),
title.size = 1,
legend.position = c(0,0),
legend.text.size=0.6,
legend.title.size = 0.6,
legend.title.fontface = 1,
legend.height = 0.4,
frame=FALSE)+
tm_credits("(d)", position=c(0,0.85), size=1)
t<-tmap_arrange(pm25_02_08_map,pm25_08_18_map,pm25_dif_map,furnaces_map, ncol=2)
t
## Variable(s) "pm25_dif" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
#tmap_save(t, 'pm25_ls.png')
